#include "Utility.h"
    
using namespace std;
void Utility::subtrBgdFromData(const TH1D *a, const TH1D *b, RooRealVar &nbgd, RooRealVar &nsig, int isol_th) {
    TH1D *all = (TH1D *) a->Clone();
    TH1D *bgd = (TH1D *) b->Clone();
    all->Print("all");
    bgd->Print("all");
    // Find bin for cut
    int bin1=0;
    for(;bin1<20;bin1++) if(floor(all->GetBinLowEdge(bin1)*100+0.5)>=isol_th) break;
    --bin1;
    double maxX = all->GetXaxis()->GetXmax();
    int    bin2 = all->FindBin(maxX);

    double int1 = all->Integral(bin1,bin2);
    double int2 = bgd->Integral(bin1,bin2);
    bgd->Scale(int1/int2);

    all->Add(bgd,-1);
    Double_t err1=-1,err2=-1;
    cout << "bin1=" << bin1 << endl;
    nbgd = bgd->IntegralAndError(0,bin2,err1);
    nbgd.setError(err1);
    nsig = all->IntegralAndError(0,bin2,err2);
    nsig.setError(err2);
    cout << "nbgd=" << nbgd.getVal() << endl;
    cout << "nsig=" << nsig.getVal() << endl;

    delete all;
    delete bgd;
}
void Utility::calcCountBelowIsolTh(const TH1D *a, const TH1D *b, RooRealVar &nbgdisolth, RooRealVar &nsigisolth, int isol_th) {
    TH1D *all = (TH1D *) a->Clone();
    TH1D *bgd = (TH1D *) b->Clone();
    int bin1=0;
    for(;bin1<20;bin1++) if(floor(all->GetBinLowEdge(bin1)*100+0.5)>=isol_th) break;
    --bin1;
    double maxX = all->GetXaxis()->GetXmax();
    int    bin2 = all->FindBin(maxX);
    double int1 = all->Integral(bin1,bin2);
    double int2 = bgd->Integral(bin1,bin2);
    bgd->Scale(int1/int2);
    all->Add(bgd,-1);
    Double_t err1=-1,err2=-1;
    nbgdisolth = bgd->IntegralAndError(0,bin1,err1);
    nbgdisolth.setError(err1);
    nsigisolth = all->IntegralAndError(0,bin1,err2);
    nsigisolth.setError(err2);
    delete all;
    delete bgd;
}


/*******************************************************************************
  Calculate integral from histogram h above/below cut and store in RooRealVar
 *******************************************************************************/
void Utility::calculateIntegral(RooRealVar &var, const TH1D* h, int cut, bool below) {
    TH1D *buf = (TH1D*) h->Clone();
    double err;

    // Find bin for cut
    int bin1=0;
    for(;bin1<20;bin1++) if(floor(buf->GetBinLowEdge(bin1)*100+0.5)>=cut) break;
    --bin1;

    int bin2 = buf->FindBin(1.5);
    if(below)
    {
        bin2=bin1;
        bin1=1;
    }

    var=buf->IntegralAndError(bin1,bin2,err);
    var.setError(err);
    delete buf;
}

/*******************************************************************************
  Correct integral to range above/below cut and store in RooRealVar
 *******************************************************************************/
void Utility::correctToRange(RooRealVar &var, TH1D* h, int cut, bool below) {
    TH1D *buf = (TH1D*) h->Clone();
    buf->Scale(var.getVal()/buf->Integral());
    double err;

    // Find bin for cut
    int bin1=0;
    for(;bin1<20;bin1++) if(floor(buf->GetBinLowEdge(bin1)*100+0.5)>=cut) break;
    --bin1;

    int bin2 = buf->FindBin(1.5);
    if(below)
    {
        bin2=bin1;
        bin1=1;
    }
    var=buf->IntegralAndError(bin1,bin2,err);
    var.setError(err); // FIX THIS ERROR!!
    delete buf;
}

